/*
===========================================================================
This source file is part of DyLab (Dynamic Laboratory)
For the latest info, see http://code.google.com/p/dylab/

Copyright (c) 2006-2008 Lukas Krejci
(krejci.lukas@volny.cz)

This file is part of DyLab.

    DyLab is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    DyLab is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with DyLab.  If not, see <http://www.gnu.org/licenses/>.
===========================================================================
*/

#ifndef __dylab_TriangleGeometry3_h__
#define __dylab_TriangleGeometry3_h__

//===========================================================================================
#include <dylab/utils/EnumIterator.h>
#include <dylab/utils/String.h>
#include <dylab/utils/geometry/Plane.h>
#include <dylab/utils/exceptions/ENumericError.h>
#include <dylab/utils/exceptions/EInvalidArgument.h>
#include <dylab/utils/exceptions/EBreak.h>

//===========================================================================================
namespace dylab {

    /** 
     */
    class DYLAB_EXPORT TriangleGeometry3 
    {
    public:
        ///
        enum vertexidx_t
        {
            V_0 = 0,
            V_1,
            V_2
        };
        ///
        class DYLAB_EXPORT VertexIdxIterator : public EnumIterator<vertexidx_t, V_0, V_2> { };
        ///
        enum { V_COUNT = V_2 + 1 };
        ///
        enum edgeidx_t
        {
            E_0 = 0,
            E_1,
            E_2
        };
        ///
        class DYLAB_EXPORT EdgeIdxIterator : public EnumIterator<edgeidx_t, E_0, E_2> { };
        ///
        enum { E_COUNT = E_2 + 1 };
        //-----------------------------------------------------------------------------
        //
        struct DYLAB_EXPORT FeatureIdx
        {
            ///
            union 
            {
                vertexidx_t vertexIdx;
                edgeidx_t   edgeIdx;
            };
            ///
            FeatureIdx()
                    { }
            ///
            FeatureIdx(vertexidx_t _vertexIdx)
                : vertexIdx(_vertexIdx)
                    { }
            ///
            FeatureIdx(edgeidx_t _edgeIdx)
                : edgeIdx(_edgeIdx)
                    { }
        };
        //-----------------------------------------------------------------------------
        class DYLAB_EXPORT pointcls_t 
            : public std::pair<Math::pointcls_t, FeatureIdx>
        { 
        public:
            pointcls_t(Math::pointcls_t pointCls, FeatureIdx featureIdx)
                : std::pair<Math::pointcls_t, FeatureIdx>(pointCls, featureIdx)
                    { }
        };
        //-----------------------------------------------------------------------------
    public:
        /**
        */
        static void getEdgeVertexIndices(edgeidx_t _edgeIdx, vertexidx_t (&_vertexIndices)[EdgeGeometry3::V_COUNT])
                {
                    switch (_edgeIdx)
                    {
                    case E_0: 
                        _vertexIndices[EdgeGeometry3::V_0] =  V_0,
                        _vertexIndices[EdgeGeometry3::V_1] =  V_1; 
                        break;
                    case E_1: 
                        _vertexIndices[EdgeGeometry3::V_0] =  V_1,
                        _vertexIndices[EdgeGeometry3::V_1] =  V_2; 
                        break;
                    case E_2: 
                        _vertexIndices[EdgeGeometry3::V_0] =  V_2,
                        _vertexIndices[EdgeGeometry3::V_1] =  V_0; 
                        break;
                    default: 
						throw EInvalidArgument(String::fmt("Invalid edge index: %1%", static_cast<uint32_t>(_edgeIdx)));
                    }
                }
    private:       
        ///
        Vector3 mVertices[V_COUNT];
    public:

        /**          
         */
        TriangleGeometry3() 
                { }

        /**
        */
        TriangleGeometry3(const Vector3 & _v0, const Vector3 & _v1, const Vector3 & _v2)             
                { mVertices[V_0] = _v0, mVertices[V_1] = _v1, mVertices[V_2] = _v2; }

        /**
        */
        TriangleGeometry3(const Vector3 (&_vertices)[V_COUNT])             
                { for (VertexIdxIterator vertexIdx; vertexIdx.isValid(); vertexIdx++) mVertices[vertexIdx] = _vertices[vertexIdx]; }

        /**
        */
        TriangleGeometry3(const TriangleGeometry3 & _triangle) 
                { for (VertexIdxIterator vertexIdx; vertexIdx.isValid(); vertexIdx++) mVertices[vertexIdx] = _triangle.mVertices[vertexIdx]; }

        /** 
         */
        Vector3 & vertex(vertexidx_t _vertexIdx) 
                { return mVertices[_vertexIdx]; }

        /** 
         */
        const Vector3 & vertex(vertexidx_t _vertexIdx) const
                { return mVertices[_vertexIdx]; }

        /**
        */
        Vector3 & operator[] (vertexidx_t _vertexIdx) 
                { return mVertices[_vertexIdx]; }
        
        /**
        */
        const Vector3 & operator[] (vertexidx_t _vertexIdx) const
                { return mVertices[_vertexIdx]; }

        /**
        */
        Vector3 & operator[] (uint32_t _index) 
                { 
                    DYLAB_ASSERT(_index < V_COUNT);
                    return mVertices[_index]; 
                }
        
        /**
        */
        const Vector3 & operator[] (uint32_t _index) const
                { 
                    DYLAB_ASSERT(_index < V_COUNT);
                    return mVertices[_index]; 
                }

        /**
        */
        Plane getPlane() const
                { return Plane(mVertices[V_0], mVertices[V_1], mVertices[V_2]); }

        /**
        */
        EdgeGeometry3 getEdgeGeometry(edgeidx_t _edgeIdx) const
                {
                    vertexidx_t edgeVertexIndices[EdgeGeometry3::V_COUNT];
                    getEdgeVertexIndices(_edgeIdx, edgeVertexIndices);
                    return EdgeGeometry3(mVertices[edgeVertexIndices[V_0]], mVertices[edgeVertexIndices[V_1]]);
                }

        /**
        */
        bool isDegenerated() const
                { return Math::almostZero(computeArea()); }

        /** 
         */
        Vector3 computeCenter() const
                { 
                    Vector3 cp = Vector3::ZERO;
                    for (VertexIdxIterator vertexIdx; vertexIdx.isValid(); vertexIdx++)
                        cp += vertex(vertexIdx);
                    return cp / V_COUNT; 
                }

        /**
        */
        Vector3 computeNormal() const
                { return (mVertices[V_1] - mVertices[V_0]).cross(mVertices[V_2] - mVertices[V_0]).normalise(); }

        /**
        */
        real_t computeArea() const
                { return (mVertices[V_1] - mVertices[V_0]).cross(mVertices[V_2] - mVertices[V_0]).length() / DYLAB_R(2); }

        /**
        */
        Vector3 computeBarycentricCoords(const Vector3 & _p) const
                {                    
                    // compute area of the triangle
                    real_t area = computeArea();

                    // shouldn't be zero
                    if (area == Math::ZERO)
                        throw ENumericError(String::fmt("Can't compute barycentric coordinates for point %1% inside the triangle %2% (degenerated: %3%)", _p, *this, isDegenerated()));

                    // compute & return the barycentric coordinates
                    Vector3 bc;
                    bc.y() = TriangleGeometry3(mVertices[V_0], mVertices[V_2], _p).computeArea() / area;
                    bc.z() = TriangleGeometry3(mVertices[V_0], mVertices[V_1], _p).computeArea() / area;
                    bc.x() = Math::ONE - bc.y() - bc.z();
                    return bc;
                }

        /**
        */
        pointcls_t classifyPoint(const Vector3 & _p) const
                {
                    // compute barycentric coordinates of the point _p inside the triangle
                    Vector3 bc = computeBarycentricCoords(_p);

                    // store the zero flag of all coords
                    bool xAlmostZero = Math::almostZero(bc.x(), DYLAB_R(100) * Math::EQUALITY_TOLERANCE);
                    bool yAlmostZero = Math::almostZero(bc.y(), DYLAB_R(100) * Math::EQUALITY_TOLERANCE);
                    bool zAlmostZero = Math::almostZero(bc.z(), DYLAB_R(100) * Math::EQUALITY_TOLERANCE);

                    // first test, whether the point is not outside
                    real_t bc2sum = bc.y() + bc.z();
                    if ((((bc.y() > Math::ZERO) || yAlmostZero) && ((bc.z() > Math::ZERO) || zAlmostZero)) && 
                          (bc2sum < Math::ONE || Math::almostEqual(bc2sum, Math::ONE)))
                    {
                        // test, if point is one of triangle's vertices
                        if (yAlmostZero && zAlmostZero)                                
                            return pointcls_t(Math::POINTCLS_VERTEX, V_0);
                        if (xAlmostZero && zAlmostZero)
                            return pointcls_t(Math::POINTCLS_VERTEX, V_1);
                        if (xAlmostZero && yAlmostZero)
                            return pointcls_t(Math::POINTCLS_VERTEX, V_2);

                        // test, if point lies in one of triangle's edges
                        if (zAlmostZero)
                            return pointcls_t(Math::POINTCLS_EDGE, E_0);
                        if (xAlmostZero)
                            return pointcls_t(Math::POINTCLS_EDGE, E_1);
                        if (yAlmostZero)
                            return pointcls_t(Math::POINTCLS_EDGE, E_2);

                        // _p is inner point
                        return pointcls_t(Math::POINTCLS_INNER, FeatureIdx());
                    }
                    else
                    {
                        // the point is outside
                        return pointcls_t(Math::POINTCLS_OUTER, FeatureIdx());
                    }
                }

        /**
        */
        Math::intersection_t computeLineIntersection(const Vector3 & _a, const Vector3 & _b, Vector3 & _ip = Vector3(), bool allowAMinus = true, bool allowBPlus = true) const
                {
                    Vector3    u, v, n;         // triangle vectors
                    Vector3    w0, w;           // ray vectors
                    real_t     r, a, b;         // params to calc ray-plane intersect                    

                    // get triangle edge vectors and plane normal
                    u = mVertices[1] - mVertices[0];
                    v = mVertices[2] - mVertices[0];
                    n = u.cross(v);  

                    // test, if triangle is degenerate
                    if (n.almostZero())         
                        return Math::INTERSECTION_NONE;  

                    // compute triangle plane intersection
                    Vector3 dir = _b - _a;
                    w0 = _a - mVertices[0];
                    a = - n.dot(w0);
                    b = n.dot(dir);

                    // test, if ray is parallel to triangle plane
                    if (Math::almostZero(b)) 
                    {      
                        // test, if ray lies in triangle plane
                        if (a == 0)                 
                            return Math::INTERSECTION_LINE;                        
                        else // ray disjoint from plane
                            return Math::INTERSECTION_NONE;           
                    }

                    // get intersect point of ray with triangle plane
                    r = a / b;

                    // test, if ray goes away from triangle, before point _a
                    if (r < Math::ZERO && !allowAMinus)        
                        return Math::INTERSECTION_NONE;                          

                    // test, ray goes after point _b
                    if (r > Math::ONE && !allowBPlus)         
                        return Math::INTERSECTION_NONE;                          

                    // intersect point of ray and plane
                    _ip = _a + dir * r;         

                    // is _ip inside Triangle?
                    real_t uu, uv, vv, wu, wv, D;
                    uu = u.dot(u);
                    uv = u.dot(v);
                    vv = v.dot(v);
                    w = _ip - mVertices[0];
                    wu = w.dot(u);
                    wv = w.dot(v);
                    D = uv * uv - uu * vv;

                    // get and test parametric coords
                    real_t s, t;

                    // test 1, if _ip is outside Triangle
                    s = (uv * wv - vv * wu) / D;
                    if (s < Math::ZERO || s > Math::ONE)       
                        return Math::INTERSECTION_NONE;
                 
                    // test 2, _ip is outside Triangle
                    t = (uv * wu - uu * wv) / D;    
                    if (t < Math::ZERO || (s + t) > Math::ONE) 
                        return Math::INTERSECTION_NONE;

                    // _ip is in Triangle
                    return Math::INTERSECTION_POINT;    
                }

        /** 
         */
        Math::intersection_t computeRayIntersection(const Ray3 & _ray, Vector3 & _ip = Vector3()) const
                { return computeLineIntersection(_ray.start(), _ray.start() + _ray.direction(), _ip, false, true); }

        /** 
         */
        Math::intersection_t computeEdgeIntersection(const EdgeGeometry3 & _edge, Vector3 & _ip = Vector3()) const
                { return computeLineIntersection(_edge.vertex(EdgeGeometry3::V_0), _edge.vertex(EdgeGeometry3::V_1), _ip, false, false); }

        /**
        */
        operator String() const                
                { return String::fmt("[%1%, %2%, %3%]", mVertices[V_0], mVertices[V_1], mVertices[V_2]); }
    };

    //===========================================================================================    
    inline std::ostream & operator << (std::ostream & _os, const TriangleGeometry3 & _o)
            { return _os << String(_o); }

    //===========================================================================================    
    inline std::ostream & operator << (std::ostream & _os, const TriangleGeometry3::pointcls_t & _pointCls)
            { 
                switch (_pointCls.first)
                {
                case Math::POINTCLS_INNER: 
                    return _os << "POINTCLS_INNER"; 
                case Math::POINTCLS_VERTEX: 
                    return _os << String::fmt("POINTCLS_VERTEX(V_%1%)", _pointCls.second.vertexIdx); 
                case Math::POINTCLS_EDGE: 
                    return _os << String::fmt("POINTCLS_EDGE(E_%1%)", _pointCls.second.edgeIdx); 
                case Math::POINTCLS_OUTER: 
                    return _os << "POINTCLS_OUTER"; 
                default: 
                    throw EInvalidArgument(String::fmt("Invalid point-triangle classification: %1%", _pointCls.first));
                    return _os;
                }
            }
}
//===========================================================================================
#endif // __dylab_TriangleGeometry3_h__
